In September 2021, a significant jump in seismic activity on the island of La Palma (Canary Islands, Spain) signaled the start of a volcanic crisis that still continues at the time of writing. Earthquake data is continually collected and published by the Instituto Geográphico Nacional (IGN). …
Rows: 110038 Columns: 3
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (3): x, y, lst_c
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
dados_latlong_sf <-st_as_sf(df, coords =c("x", "y"), crs =4326)# 2. Reprojetar para o mesmo CRS do municípiodados_latlong_sf_proj <-st_transform(dados_latlong_sf, crs =st_crs(ride_sf))# 3. Filtrar apenas os pontos *dentro* do municípiodados_ride_latlong <-st_filter(dados_latlong_sf_proj, ride_sf) # ✅ correto!dados_ride_latlong <- dados_ride_latlong[seq(1, nrow(dados_ride_latlong), by =18), ] #15 x18 16 # 6. Criar banco com latitude/longitude (graus decimais)dados_latlong <- dados_ride_latlong %>%mutate(x =st_coordinates(geometry)[,1],y =st_coordinates(geometry)[,2] ) %>%select(lst_c, x, y) %>%as.data.frame() %>%select(-geometry) %>%rename(avg_z=lst_c)# 7. Transformar para UTM (zona 23S - SIRGAS 2000 / EPSG:31983)dados_utm_sf <-st_transform(dados_ride_latlong, crs =31983)# 8. Criar banco com coordenadas UTM (em metros)dados_utm <- dados_utm_sf %>%mutate(x =st_coordinates(geometry)[,1],y =st_coordinates(geometry)[,2] ) %>%select(lst_c, x, y) %>%as.data.frame() %>%select(-geometry)%>%rename(avg_z=lst_c)# 9. Exportar os dois arquivoswrite.csv(dados_latlong, "RIDE_temp_latlong.csv", row.names =FALSE)write.csv(dados_utm, "RIDE_temp_utm.csv", row.names =FALSE)#set.seed(181326) #set.seed(181019)# para reprodutibilidadeset.seed(020279)dados_utm_sf_amostra <- dados_utm_sf |>rename(temp = lst_c) |>sample_n(150)
In [6]:
ggplot() +geom_sf(data = ride_sf, fill =NA, color ="black", linewidth =0.1) +# polígono em pretogeom_sf(data = dados_utm_sf_amostra, aes(color = temp), size =2) +# pontos coloridosscale_color_viridis_c(name ="Temperatura (°C)") +theme_minimal() +labs(title ="Temperatura na região da RIDE",subtitle ="Valores extraídos de MODIS (LST)",caption ="Fonte: MODIS via modisfast" )
In [7]:
v <-variogram(temp ~1, data = dados_utm_sf_amostra)v$dist <- v$dist /1000# converte distância para kmfv <-fit.variogram(v, model =vgm(psill =15, model ="Sph", range =200, nugget =3))plot(v, fv)
fv <-fit.variogram(v, model =vgm(psill =15, model ="Gau", range =200, nugget =3))plot(v, fv)
fv <-fit.variogram(v, model =vgm(psill =10, model ="Exp", range =200, nugget =2))plot(v, fv)
In [8]:
# Transformar ambos para UTM zona 23Sride_utm <-st_transform(ride_sf, 31983)# ✅ Verificação da variável de temperaturanames(dados_utm_sf_amostra)
[1] "temp" "geometry"
summary(dados_utm_sf_amostra$temp)
Min. 1st Qu. Median Mean 3rd Qu. Max.
27.04 31.72 33.41 33.99 35.65 44.91
# 📈 Variograma empírico e ajustevc <-variogram(temp ~1, dados_utm_sf_amostra)vinitial <-vgm(psill =15, model ="Exp", range =150000, nugget =3)fv <-fit.variogram(vc, model = vinitial)plot(vc, fv)
# 🔳 Criar grade regular de pontos dentro da RIDEgrid <-st_make_grid(ride_utm, cellsize =1000, what ="centers")grid_sf <-st_sf(geometry = grid, crs =st_crs(ride_utm))grid_sf <- grid_sf[ride_utm, ]# 🤖 Krigagemk <-gstat(formula = temp ~1, data = dados_utm_sf_amostra, model = fv)kpred <-predict(k, newdata = grid_sf)
[using ordinary kriging]
kpred_sf <-st_as_sf(kpred)# 🌍 Projeção para latitude/longituderide_latlon <-st_transform(ride_utm, 4326)kpred_latlon <-st_transform(kpred_sf, 4326)# 🎨 Plot final com paleta magma + nomes dos municípios\ggplot() +geom_sf(data = kpred_latlon, aes(color = var1.pred), size =0.6) +geom_sf(data = ride_latlon, fill =NA, color ="black",linewidth =0.1) +geom_sf_text(data = ride_latlon, aes(label = NM_MUNICIP), size =2, color ="black")+# scale_color_viridis_c(name ="Temperatura prevista (°C)", option ="magma") +theme_minimal() +labs(title ="Krigagem da Temperatura MODIS",subtitle ="Interpolação espacial para a RIDE-DF",caption ="Fonte: MODIS via modisfast" )
Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
give correct results for longitude/latitude data
In [9]:
# 📦 Pacotesif (!require("pacman")) install.packages("pacman")pacman::p_load(sf, gstat, stars, leaflet, viridis)# 🌍 Reprojetar ambos para WGS84 (latitude/longitude)ride_ll <-st_transform(ride_utm, 4326)kpred_ll <-st_transform(st_as_sf(kpred), 4326)# 🎨 Paleta de cores (magma)pal <-colorNumeric(palette ="magma", domain = kpred_ll$var1.pred)# 🗺️ Mapa interativo com limites e nomes dos municípiosleaflet() %>%addProviderTiles(providers$CartoDB.Positron) %>%addCircleMarkers(data = kpred_ll,radius =3,fillColor =~pal(var1.pred),fillOpacity =0.7,stroke =FALSE,label =~paste0("Temp: ", round(var1.pred, 1), " °C") ) %>%# Polígono com os municípios (com rótulo interativo)addPolygons(data = ride_ll,color ="white",weight =1,fillOpacity =0,label =~NM_MUNICIP # altere aqui caso o nome esteja em outra coluna ) %>%# Limite da RIDE destacado em pretoaddPolygons(data =st_union(ride_ll), # união dos polígonos para contorno da RIDEcolor ="black",weight =2,fill =FALSE,label ="Limite da RIDE" ) %>%# Pontos de krigagem coloridos# LegendaaddLegend(position ="bottomright",pal = pal,values = kpred_ll$var1.pred,title ="Temperatura prevista (°C)",opacity =1 )